Phase transition classes in triplet and quadruplet reaction diffusion models 
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Phase transitions of reaction-diffusion systems with site occupation restriction and with particle 
creation that requires n = 3, 4 parents, whereas explicit diffusion of single particles (A) is present 
are investigated in low dimensions by mean-field approximation and simulations. The mean-field 
approximation of general nA — > (n + k)A, mA — > (rn — I) A type of lattice models is solved and novel 
kind of critical behavior is pointed out. In d = 2 dimensions the 3^4 — > AA, 3A — > 2A model exhibits 
a continuous mean-field type of phase transition, that implies d c < 2 upper critical dimension. 
For this model in d = 1 extensive simulations support a mean-field type of phase transition with 
logarithmic corrections unlike the Park et al.'s recent study (Phys. Rev E 66, 025I0I (2002)). On 
the other hand the AA — > 5A, 4A — » 3A quadruplet model exhibits a mean-field type of phase 
transition with logarithmic corrections in d = 2, while quadruplet models in Id show robust, non- 
trivial transitions suggesting d c — 2. Furthermore I show that a parity conserving model 3A —* 5A, 
2A — > in d — 1 has a continuous phase transition with novel kind of exponents. These results are 
in contradiction with the recently suggested implications of a phenomenological, multiplicative noise 
Langevin equation approach and with the simulations on suppressed bosonic systems by Kockelkoren 
and Chate (cond- mat/0208497). 



I. INTRODUCTION 

Phase transitions in genuine nonequilibrium systems 
have been investigated often among reaction-diffusion 
(RD) type of models exhibiting absorbing states [1-3]. 
In many cases mapping to surface growth, spin systems 
or stochastic cellular automata can be done. The classi- 
fication of universality classes of second order transitions 
is still one of the most important uncompleted task. One 
hopes that symmetries and spatial dimensions are the 
most significant ingredients as in equilibrium cases, how- 
ever it turned out that in many cases there is a short- 
age of such factors to explain novel universality classes. 
An important example was being investigated during the 
past two years that emerges at phase transitions of bi- 
nary production systems [4-13] (PCPD). In these sys- 
tems particle production competes with pair annihila- 
tion and single particle diffusion. If the production wins 
steady states with finite particle density appear in (site 
restricted) models with hard-core repulsion, while in un- 
restricted (bosonic) models the density diverges. By low- 
ering the production/annihilation rate a doublet of ab- 
sorbing states without symmetries emerges. One of such 
states is completely empty, the other possesses a single 
wandering particle. In case of site restricted systems the 
transition to absorbing states is continuous. 

Although the nature of this transition has not com- 
pletely been settled numerically and by field theory yet, 
an other novel class appearing in triplet production sys- 
tems was proposed very recently [14,15] (TCPD). This 
reaction-diffusion model differs from the PCPD that for 
new particle generation at least three particles have to 
meet. It is important to note that these models do not 
break the DP hypothesis [16,17] — - according to which 



in one component systems exhibiting continuous phase 
transitions to single absorbing state (without extra sym- 
metry and inhomogeneity or disorder) short ranged in- 
teractions can generate DP class transition only — - be- 
cause they exhibit multiple absorbing states that are not 
frozen, lonely particle (s) may diffuse in them. 

A phenomenologically introduced Langevin equation 
that exhibits real, multiplicative noise was suggested 
[14] to describe the critical behavior of reaction-diffusion 
models of types 



nA — > (n + 1)^, nA — > jA, 



(1) 



(with j < n number of interacting particles) in the form 

d t p{x, t) = ap(x, t) n - p(x, t) n+1 + DV 2 p(x, t) + ((x, t), 

(2) 



with noise correlations 



< ({x, t)((x',t') >= Tp fl S d (x - x')8{t - t'). 



(3) 



The classification of universality classes of nonequilib- 
rium systems by the exponent p of a multiplicative noise 
in the Langevin equation was suggested some time ago 
by Grinstein et al. [18]. However it turned out that there 
may not be corresponding particle systems to real mul- 
tiplicative noise cases [4] and an imaginary part appears 
as well if one derives the Langevin equation of a RD sys- 
tem starting from the Master equation in a proper way. 
This observation led Howard and Tauber to investigate 
systems with complex noise appearing in binary produc- 
tion models. Unfortunately the cases with and without 
occupation number restriction turned out to be different 
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in d = 1, although in d = 2 this difference was found to 
disappear at criticality and below [13]. 

By rescaling eq.(2) one can get the corresponding 
mean-field critical exponents 



(3 Mt = 1, 



MF 



= n/2, 



u!f F = n. 



(4) 



The authors of [14] expect that the noise exponent should 
be in the range 



1 < p < n, 



(5) 



hence by simple power counting the upper critical dimen- 
sion should be 



d c = 2 + 



4-2p 



(6) 



This implies for a triplet processes: 4/3 < d c < 8/3 and 
for a quadruplet (n = 4) processes: 1 < d c < 5/2. 

Very recently Kockelkoren and Chatc introduced 
stochastic cellular automata (SCA) versions of general 
nA — > (n + k)A, mA — > (m — I) A type of models [15], 
where multiple particle creation on a given site is sup- 
pressed by an exponentially decreasing creation probabil- 
ity (p N / 2 ) of the particle number. They claim that their 
simulation results in Id are in agreement with the fully 
occupation number restriction counterparts and set up a 
general table of universality classes, where as the func- 
tion of n and m only 4 classes exist, namely the directed 
percolation class [16,17], the parity conserving class [19], 
the PCPD and TCPD classes. 

In any case the heuristic Langevin equation with real 
noise assumption for RD models [14,15] should be proven 
for n > 1. Furthermore in low dimensions topological 
constraints may cause different critical behavior with and 
without occupation number restriction [20] . Note that in 
case of binary production models it had not been clear at 
all if the d c = 2 prediction of the bosonic field theory had 
also been true for site restricted systems until the numer- 
ical confirmation of [13]. In this paper I show simulation 
results for lattice models with restricted site occupancy 
in d = 1,2 with the aim of locating the upper critical 
dimensions and checking claims the of refs. [14,15] about 
possible new universality classes. 



II. MEAN-FIELD CONSIDERATIONS 

In this section I discuss the mean-field equation that 
can be set up for site restricted lattice models with gen- 
eral microscopic processes of the form 



iA^(n + k)A, mA±(m-l)A, 



(7) 



with n > 1, m > 1, k > 0, I > and m - I > 0. Note 
that this formulation is different from that of eq.(2) that 
is valid for coarse grained, continuous bosonic description 
of these reaction-diffusion systems. In this case the dif- 
fusion drops out and one can neglect the noise, hence the 



the competition of creation (with probability oa) and an- 
nihilation or coagulation (parametrized with probability 
A = 1 — a) is left behind 



(8) 



where p denotes the site occupancy probability and a is a 
dimension dependent coordination number. Each empty 
site has a probability (l-p) in mean-field approximation, 
hence the need for k empty sites at a creation brings in 
a (1 — p) k probability factor. By expanding (1 — p) k and 
keeping the lowest order contribution one can see that 
for site restricted lattice systems a p n+1 -th order term 
appears automatically with negative coefficient that reg- 
ulates eq.(8). The steady state solution can be found 
analytically in many cases and may result in different, 
continuous or discontinuous phase transitions. Here I 
split the discussion of the solutions to three parts: (a) 
n = m, (b) n > m and (c) n < m. In the inactive 
phases one expects a dynamical behavior described by 
the mA — > process, for which p oc t 1 /( m ~ 1 ) is known 
[19]. 



A. The n — m symmetric case 

The steady state solution in this case can be obtained 
by solving 



ku{\- pf = l{l-u), 



(9) 



where the trivial (p = 0) solution has been factored out. 
For the active phase one gets 



p=l- 

which vanishes at a c = 
larity 



ll-a 
k a 



l/k 



(10) 



k+i 



with the leading order singu- 



p oc |<7 - <7 C 



(11) 



and order parameter exponent exponent f3 MF = 1. At 
the critical point the time dependent behavior is de- 
scribed by 



^ = -2afcV +1 + 0(p"+ 2 ). 
that gives a leading order power-law solution 



(12) 



p cx r x i n (13) 

ned 

bosonic, coarse grained formulation in [14] too. 



hence auF — (i MF jv^ F = l/n. This was obtained from 
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B. The n > m case 



In this case besides the p = absorbing state solution 
we can get an active state if 



kap n ~ m {\- p) k =1(1- a) 



(14) 



is satisfied. Both sides are linear functions of a such that 
for (7 — ► only the p = is a solution. The left hand 
side is a convex function of p (from above) with zeros at 
p = and p = 1 . 




FIG. 1. Steady state mean-field solution for (a) n > m and 

(b) n < m cases 



Therefore by increasing a from zero the left hand side 
meets the right hand side at a c , p c > (See Fig. 1(a)). If 
this solution is stable a first order transition takes place in 
the system. Note that in higher order cluster mean-field 
solutions, where the diffusion can play a role the transi- 
tion may turn into continuous one [22-24] , therefore it is 
important to check the type of transition for d > d c . In 
Section III C I shall confirm the first orderedness of such 
transitions for two models in 2d. 



C. The n < m case 
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FIG. 2. Density times t i/3 in the two dimensional 
3A -» 4A, 3A -» 2A model for D = 0.5 and p = 0.496, 
0.4965, 0.4967, 0.4968, 0.497, 0.498, 0.4985 0.499 (top to bot- 
tom curves). The insert shows the corresponding local slopes. 



III. SIMULATIONS IN TWO DIMENSIONS 

Two dimensional simulations were performed on L = 
400 — 1000 linear sized lattices with periodic boundary 
conditions. One Monte Carlo step (MCS) — correspond- 
ing to dt — l/P (where P is the number of particles) - 
is built up from the following processes. A particle and 
a random number x G (0,1) are selected randomly; if 
x < D a site exchange is attempted with one of the ran- 
domly selected empty nearest neighbors (nn); if x > D 
k number of new particles are created with probability 
(1 — p) at randomly selected empty nn sites provided the 
number of nn particles was greater than or equal n; or 
if x > D I number of particles are removed with proba- 
bility p (taking into account the m — I > condition as 
well). The simulations were started from fully occupied 
lattices and the particle density decay was measured up 
to 10 6 - 10 8 MCS. 



By factoring out the trivial p = solution we are faced 
with the general condition for a steady state 



ka(l - p) k = 1(1 - a)p n 



(15) 



One can easily check that in this case the critical point is 
at a c = (see Fig. 1(b)) and here the density decays with 
ct = l/(m — 1) as in case of the n = 1 branching and 
to = / annihilating models showed by Cardy and Tauber 
[19] (BkARW classes). However the steady state solu- 
tion for ii > 1 gives different (3 exponents than those of 
BkARW classes, namely (3 MF = l/(m — n). This imply 
novel kind of critical behavior in low dimensions, that 
should be a subject of further investigations [21]. 



A. The 3A -> 4A, 3A — > 2A symmetric triplet model 

First I checked the dynamic behavior in the inactive 
phase for D = 0.5 diffusion rate. At p = 0.9 one can see 
the appearance of the mean-field behavior p(t) oc i -1 / 2 
following 2 x 10 6 MCS. By decreasing p this scaling sets 
in later and later times. As Fig. 2 shows for L — 1000 sys- 
tems with t max = 10 7 MCS curves with p < 0.4965 veer 
up - corresponding to the active phase — while curves 
with p > 0.497 veer down - corresponding to the absorb- 
ing state. From the p(t) data I determined the effective 
exponents (the local slopes) defined as 



(*) = 



\n[p(t)/p(t/m)} 
ln(TO) 



(16) 
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(where I used m = 4). The critical point is estimated 
at p = 0.4967(2) with a = 0.33(1) (for local slopes see 
insert of Fig. 2). This value agrees with the mean-field 
value a MF = 1/3. 
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FIG. 3. Effective order parameter exponent results in 2d. 
Bullets correspond to the 3A — » 4A, 3A — > 2A model; squares 
to the 4A -» 5A, 4A -» 3A model at D = 0.5. The insert 
shows the logarithmic fitting for the 4A — > 5A, 4A — > 3A 
model. 
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FIG. 4. Density times t 1 ' 4 in the two dimensional 
AA -> 5A, 4A -» model for D = 0.5 and p = 0.469, 
0.47, 0.4792, 0.4705, 0.471, 0.4715, 0.4725, 0.473, 0.474, 0.476, 
0.478, 0.48 (top to bottom curves). The insert shows the cor- 
responding local slopes. 



B. The 4A — ► 5 A, 4A — » 3 A symmetric quadruplet 
model 



Density decays for several p-s in the active phase 
(0.003 < e = \p c — p\ < 0.3) were followed on logarithmic 
time scales and averaging was done over ~ 100 indepen- 
dent runs in a time window, which exceeds the level-off 
time by a decade. The steady state density in the active 
phase at a critical phase transition is expected to scale 
as 



p(oo,p) oc \p-p c f 



(17) 



Using the local slopes method one can get a precise esti- 
mate for (3 as well as for the corrections to scaling 



lnp(oo,e. t ) - lnp(oo,ej_i) 
ln(ei) - ln(ei_i) 



(18) 



Here simulations are much slower than in case of the 
triplet model, hence systems with linear size L = 400 
could be investigated. First I checked the dynamic be- 
havior in the inactive phase for D = 0.5. At p = 0.9 
a mean-field type of decay p(t) oc i" 1 / 3 can be ob- 
served following 10 6 MCS. As one can see on Fig. 4 for 
p < 0.4702 the density decay curves veer up, while for 
p > 0.4705 they veer down. The estimated critical point 
is p c ~ 0.4703(1). The effective exponent at p c extrapo- 
lates to a ~ 0.215(5). As one can see on this graph the 
separatrix (critical) curve exhibits a linear shape on the 
p(i)i 1//4 - ln(t) scale suggesting logarithmic corrections 
to scaling. Similarly the effective exponents of (3 seem to 
extrapolate to j3 ~ 0.71(5) (Fig. 3) that is very far from 
the mean- field value p MF = 1. To check the possibil- 
ity that a logarithmic correction can result in mean-field 
exponents the fitting with the lowest order correction 



p{oo,p) = [(p c -p)(a + bln(p c -p))} 
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(19) 



where I used the p c value determined before. One can 
see on Fig. 3 that the effective exponent for e > 0.005 
exhibits a correction to scaling (inclined line) and tends 
to lim e ^o P — 1-0(1), which agrees with the mean-field 
value again. By neither the a, nor the (3 exponent can 
one observe logarithmic corrections suggesting d c < 2. 

The density decay simulations were repeated at D = 
0.2, where the critical point was found at p c = 0.4795(1) 
with mean-field like a exponent again. 



has been applied for the steady state p(oo,p) data. I 
used the non-linear fitting of the "xmgr" graphical pack- 
age with a relative error in the sum of squares with at 
most 0.0001. This resulted in (3 = 1.01 at p = 0.471 
(a = -10.8, b = -6.05) (see insert of Fig. 3). This result 
in agreement with the dynamical scaling conclusion may 
support that the upper critical dimension for quadruplet 
models is rf c = 2. To get more solid results further, very 
extensive simulations would be necessary that is beyond 
the scope of this study. In any case clear mean-field be- 
havior can't be concluded. 
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FIG. 5. Density decay in the two dimensional AA — > 5A, 
2A^% model for D = 0.5 and p = 0.05, 0.119, 0.121, 0.122, 
0.123, 0.124, 0.125, 0.13 (top to bottom curves) with system 
sizes L = 400. 



C. 3A -> 5A, 2A -» and AA -> 5A, 2A -» hybrid 
models 

One can find two regions in the density decay behav- 
ior by varying p in both models. For p < p c steady 
state values are reached quickly while for p > p c a rapid 
(faster than power-law) initial density decay crosses over 
to p oc t~ x . This is in agreement with the mean-field 
behavior of the 2A — > process in Id [25] dominating in 
the inactive phase. For the AA — * 5A, 2A — » quadru- 
plet production model this threshold is at p c = 0.119(1) 
(see Fig. 5) where an abrupt jump is observable from 
p{oo) = 0.833 to p{oo) = 0. In case of the 3A — > 5A, 
2A — > triplet production model the threshold is at 
p c = 0.220(1) with a jump from p(oo) = 0.45 to zero. 

In neither cases do we see dynamical scaling at the 
transition. These results are in agreement with the first 
order transition of the mean-field approximations given 
in Sect. II B for n > m. 



IV. SIMULATIONS IN ONE DIMENSION 

The simulations in one dimension were carried out on 
L = 10 5 sized systems with periodic boundary condi- 
tions. The initial states were again fully occupied lat- 
tices, and the density of particles is followed up to 10 6 
MCS. An elementary MCS consists of the following pro- 
cesses: 

(a) A% <-> A with probability D, 

(b) mA — > (m — I) A with probability p(l — D), 

(c) nA — > (n + k)A with probability (1 — p)(l — D), 

such that the reactions were allowed on the left or right 
side of the selected particle strings. 



A. 3A -> AA, 3A 2A and 3A -> 6A, 3A -> symmetric 
triplet models 



The 3A — > AA, 3A -> 2A site restricted model in 
Id was simulated by Park et al. [14] for small sys- 
tems up to 10 6 MCS. They concluded to find a novel 
kind of phase transition with the order parameter ex- 
ponents a = 0.32(1) and (3 = 0.78(3). For the restricted 
bosonic version of this model large scale simulations gave 
a = 0.27(1) and = 0.90(5) [15]. Note that since reac- 
tive and diffusive sectors arise in this model like in PCPD, 
diffusion dependence or corrections to scaling may ham- 
per to see real asymptotic behavior [7,39,40]. Here I 
show extended simulation results for the strictly site re- 
stricted lattice model model with t max = 10 7 MCS at 
diffusion rate D — 0.1. At the critical point the a e ff(t) 
curve exhibits a straight line shape for t — > oc, while 
in sub(super)-critical cases a e //(t) curves veer down(up) 
respectively. As one can see on Fig. 6 following a long 
relaxation p < 0.3032 curves veer up, while p > 0.3035 
curves veer down in the t — > oo limit. From this one 
can estimate p — p c ~ 0.3033(1) with a = 0.33(1) in 
agreement with the results of [14]. 

By analyzing super-critical, steady state densities with 
the local slopes method one can read-off: (3 e ff — + /3 ~ 
0.95(5) (see Fig. 7), which is higher than the results of 
[14] and [15]. 

0.18 | 1 1 1 1 1 1 1 1 1 -> 




FIG. 6. Local slopes of the density decay in the Id 
3 A -> A A, 3 A -» A A model at D = 0.1. Different curves 
correspond to p = 0.3, 0.301, 0.3015, 0.302 0.3025, 0.3027, 
0.303, 0.3035, 0.304 and 0.3045 (from top to bottom). 

However one should be careful and check diffusion de- 
pendence and corrections to scaling especially because 
these critical exponent estimates are quite close to the 
mean-field values (a MF = 1/3, (3 MF = 1) and as it was 
shown in Sect. Ill d c < 2. Since the p = 0.303 and 
p = 0.3035 p(t) curves show clear curvature for large 
times the 0.303 < p c < 0.3035 conditions seems to be 
inevitable. 
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FIG. 7. Effective order parameter exponent results in Id. 
Stars correspond to 3 A — > 44, 3 A — > 24 model; bullets to 
44 -» 54, 24 -» model at 73 = 0.2; squares to 44 -» 5,4, 
24 -» model at D = 0.8; diamonds to 44 -> 54, 44 -> 
model at D = 0.3. 



[10,13]. Therefore it is still an open question whether 
parity conservation is relevant in other models than in 
BARW types. 
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FIG. 8. Spatio-temporal evolution of the critical, 1+1 d 
diffusive 3 A -> 54, 24 ^ model (73 = 0.8). 



I tried to fit the steady state data in the 0.303 < p < 
0.3035 region by the logarithmic correction form (19) and 
obtained /3 = 1.07(10) at p c — 0.3032 that agrees with 
the mean-field value and implies d c = 1. 

Just considering mean-field results, according to which 
k does not play a role for n = m models one may ex- 
pect that the 34 — > 64, 34 — > triplet creation model 
exhibits the same kind of transition as the 34 — > 44, 
34 — > 24 model. Indeed Kockelkoren and Chafe's simu- 
lations show this [15]. However doing lattice simulations 
with site restrictions it turned out that the 34 — > 64 cre- 
ation was so effective that it shifted the transition to the 
zero production limit (p — 1) where the 34 — > process 
in Id is known to decay as p oc (ln(£) /t) 2 [19] . Off-critical 
simulations showed that j3 — 0.33(1), meaning that this 
transition belongs to the BkARW mean-field class. On 
the other hand there may be other realizations of this 
model, where the transition reported by Kockelkoren and 
Chate is accessible. 



B. The 34 -> 54, 24 -» model 

It has been established that in n = 1, m = I = 2 and 
even k - so called even number of offspringed branch- 
ing and annihilating models (BARWe) - the parity con- 
serving (PC) class continuous phase transition emerges 
[19,27,28]. This class has also been observed in mod- 
els exhibiting Z2 symmetric absorbing states, where the 
domain walls separating ordered phases follow BARWe 
dynamics [29-32,35]. This class was originally called 
parity conserving, owing to the conservation law that 
made it different from the robust DP class. However 
it turned out that in Z2 symmetric models this conserva- 
tion is not enough [33-36] . Furthermore in binary spread- 
ing models this conservation was found to be irrelevant 



I investigated the phase transition of the triplet pro- 
duction 34 —¥ 54, 24 — > model (with explicit particle 
diffusion) possessing parity conservation. As I showed in 
Sect. Ill C in two dimensions this system exhibits a first 
order transition in agreement with the mean-field results. 
This first order mean-field behavior does not give a di- 
rect hint on the type of phase transition in Id. Kockelko- 
ren and Chafe's simulations on the one dimensional, sup- 
pressed bosonic cellular automaton version of this model 
shows simple DP class density decay [15]. However if we 
consider the space-time evolution we see very non-DP like 
spatio-temporal pattern (see Fig. 8). This pattern resem- 
bles much more to those of the PCPD class models, where 
compact domains separated by clouds of lonely wander- 
ing particles occur. Of course such qualitative judgment 
on the universal behavior is not enough but has been 
found to be quite successful in case of binary production 
systems [8,11]. 

The density decay simulations at D — 0.8 and D = 
0.2 have been analyzed by the local slopes method see 
Fig. 9. At D = 0.8 the critical point is estimated at 
pH = 0.4629(3) and the corresponding effective expo- 
nent tends to a H = 0.24(1). At D = 0.2-s the critical 
point is at p^ = 0.2240(3), and the local exponent seems 
to extrapolate to a L = 0.28(1). Such small difference 
between the high and low D results can also be observed 
by analyzing the steady state results: f3 H = 0.43(3) ver- 
sus P L = 0.63(3). These exponent estimates are far from 
the 1+1 d DP values (a = 0.159464(6), (3 = 0.276486(8) 
[42]), hence the claim of Kockelkoren and Chate for the 
critical behavior of n > m models is questionable. 








„ 0.23 -,Vr 



FIG. 9. Local slopes of the density decay in the Id 
3A -> 5A, 2A -» model for D = 0.8. Different curves 
correspond to p = 0.46, 0.461, 0.462, 0.4625, 0.4627, 0.463, 
0.46325, 0.464 and 0.465 (from top to bottom). 

On the other hand the diffusion dependence of the crit- 
ical exponents is a challenge and has been observed in 
the binary production PCPD model [7]. In [40] it was 
shown that assuming logarithmic corrections to scaling - 
that is quite common in Id models - a single universality 
class can be supported numerically. Therefore here again 
I have investigated the possibility of the collapse of the 
high and low D exponents. Assuming the same kind of 
logarithmic correction forms as in [40] 



[a + b]n(t)]/t] a 



(20) 



I have found a consistent set of exponents both for 
D = 0.2 and D = 0.8 (see Table I and insert of Fig.9). 
For the data analysis I used non-linear fitting of the pro- 
gram "xmgr" package, with a relative error in the sum of 
squares with at most 0.001. 



a = 0.22(1), = 0.60(1) 



(21) 



with the critical thresholds: pf = 0.4627(1), p L c = 
0.2240(1). These exponents suggest a distinct univer- 
sality class from the known ones [3] . 



C. AA -> 5A, AA 



and 4A -> 5A, 2A 
models 



quadruplet 



Two dimensional simulations (Sect. Ill), showed that 
for n = m = 4 symmetric quadruplet models d c = 
2. Simulations in the corresponding suppressed bosonic 
SCA [15] with n = 4 and 1 < m < 4 located the 
phase transition at zero production rate. Here I show 
that in the one dimensional 4 A — » 5 A, 4 A — > and 
4A — > 5A, 2A — * site restricted models continuous 
phase transitions with p < 1 and with non-trivial expo- 
nents can be found. The density decay was followed up 
to t = 10 6 MCS and the critical point was located by 
the local-slopes method (see Fig. 10) at p = 0.9028(1) for 
D = 0.3. The corresponding exponent can be estimated 



as a = 0.27(1). For D = 0.05 one gets p c = 0.9605(3) 
with a = 0.28(1), so one can not see diffusion depen- 
dence here. Analyzing off-critical data with the local 
slopes method (18) one gets (3 = 0.48(2) (see Fig.7). 



0.23 
0.24 



0.25 
0.26 




FIG. 10. Local slopes of the density decay in the Id 
4A — > 5,4, 4A — » model for D = 0.3. Different curves 
correspond to p = 0.904, 0.9037, 0.9033, 0.903 0.9027, 0.902 
(from bottom to top). 

In accordance with these results simulations for the 
4A -> 5A, 2A -> model at D = 0.2 and D = 0.8 dif- 
fusion rates resulted in p c (0.2) = 0.53185(5), p c (0.8) = 
0.5742(1) with a = 0.27(1) and (3 = 0.48(2) exponents 
(see Fig.7). As we can see critical exponent data for 
quadruplet models are robust and no diffusion depen- 
dence has been found. Furthermore critical space-time 
plots are very similar to that of the PCPD model. 



V. CONCLUSIONS 

In this paper I investigated the phase transitions of 
general nA — > (n + k)A, mA — > (m — I) A reaction type 
of models with explicit single particle diffusion on oc- 
cupation number restricted lattices in one and two di- 
mensions. I showed that mean-field solution for n — m 
symmetric cases results in universality classes character- 
ized by the exponents a = 1/n, fi = 1. I determined the 
upper critical dimensions for the triplet and quadruplet 
cases by simulations. For n — 3 high precision simula- 
tions show mean-field type of criticality with logarithmic 
corrections meaning d c — 1. This result is in contradic- 
tion with the simulations of [14] and [15] and with the 
analytical form for d c (n) derived from a phenomenologi- 
cal Langevin equation. In case of my site restricted real- 
ization of the one dimensional 3A — > 6A, 3A — > model 
the phase transition point is shifted to zero production 
rate and is continuous, BkARW mean-field type. This 
is in contradiction with the findings of [15] for an other 
stochastic cellular automaton realization of this model. 
For n — 4 the upper critical dimension was located at 
d c = 2 opening up the possibility for non-trivial critical 
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behavior in d = 1. Indeed two versions of such quadru- 
plet models were shown to exhibit robust, novel type of 
critical transition in one dimension. 

For n > m the mean-field approximation gives first or- 
der transition that was observed by simulations for two 
(n = 3,4) models in d = 2. On the other hand numeri- 
cal evidence was given that the parity conserving model 
3A — > 5A, 2A — > in Id exhibits a non-PC type of criti- 
cal behavior with logarithmic corrections by varying the 
diffusion rate. This transition does not fit in the univer- 
sality class scheme suggested by [15]. 

Finally I showed that for n < m models the mean- 
field approximations result in new classes featured by 
a MF = l/(m - 1) and (3 MF = l/(m - n). Such kind 
of models should be subject of further studies. 

The presented mean-field and simulation results show 
that the universal behavior of such low-dimensional 
reaction-diffusion models is rich and the table of univer- 
sality classes given by [15] is not valid for Id, fully site re- 
stricted systems. Perhaps the strict site restriction plays 
an important role that causes the differences. Field the- 
oretical (possibly fermionic) treatment that starts from 
the master equation should be set up to determine at 
least the analytical form of d c (n) for n = m > 2 models. 

Acknowledgements: 

I thank H. Chate, H. Hinrichsen and U. Tauber for use- 
ful communications and N. Menyhard for her comments 
to the manuscript. Support from Hungarian research 
funds OTKA (Grant No. T-25286), Bolyai (Grant No. 
BO/00142/99) and IKTA (Project No. 00111/2000) is 
acknowledged. The simulations were performed on the 
parallel cluster of SZTAKI and on the supercomputer 
of NIIF Hungary within the scope of the DEMOGRID 
project. 



[1] For references see : J. Marro and R. Dickman, Nonequi- 
hbrium phase transitions in lattice models, Cambridge 
University Press, Cambridge, 1999. 

[2] For a review see : H. Hinrichsen, Adv. Phys. 49, 815 
(2000). 

[3] For a recent review see : G. Odor, Phase transition 
universality classes of classical, nonequilibrium systems, 
cond-mat/0205644. 

[4] M. J. Howard and U. C. Tauber, J. Phys. A 30, 7721 
(1997). 

[5] E. Carlon, M. Henkel and U. Schollwock, Phys. Rev. E 

63, 036101-1 (2001). 
[6] H. Hinrichsen, Phys. Rev. E 63, 036102-1 (2001). 
[7] G. Odor , Phys. Rev. E62, R3027 (2000). 
[8] H. Hinrichsen, Physica A 291, 275-286 (2001). 



[9] G. Odor , Phys. Rev. E 63, 067104 (2001). 
[10] K. Park, H. Hinrichsen, and In-mook Kim, Phys. Rev. E 

63, 065103(R) (2001). 
[11] G. Odor, Phys. Rev. E 65, 026121 (2002). 
[12] J. D. Noh and H. Park, cond-mat/0109516. 
[13] G. Odor, M. A. Santos, M. C. Marques, Phys. Rev. E 

65, 056113 (2002). 
[14] K. Park, H. Hinrichsen and I. Kim, Phys. Rev. E 66, 

025101 (2002). 
[15] J. Kockelkoren and H. Chate, cond-mat/0208497. 
[16] H. K. Janssen, Z. Phys. B 42, 151 (1981). 
[17] P. Grassberger, Z. Phys. B 47, 365 (1982). 
[18] Grinstein G, Munoz M.A., and Tu Y., Phys. Rev. Lett. 

76, 4376 (1996); Tu Y., Grinstein G, and Munoz M.A. 

Phys. Rev. Lett. 78, 274 (1997). 
[19] J. L. Cardy and U. C. Tauber, Phys. Rev. Lett. 77, 4780 

(1996); J. L. Cardy and U. C. Tauber, J. Stat. Phys. 90, 

1 (1998). 

[20] G. Odor and N. Menyhard, Physica D 168, 305 (2002). 
[21] G. Odor, cond-mat/0304023. 

[22] G. Odor, N. Boccara, G. Szabo, Phys. Rev. E 48, 3168 

(1993) . 

[23] G. Odor and A. Szolnoki, Phys. Rev. E 53, 2231 (1996). 
[24] N. Menyhard and G. Odor, J. Phys. A. 28, 4505 (1995). 
[25] B. P. Lee, J. Phys. A 27, 2633 (1994). 
[26] F. van Wijland, K. Oerding and H. J. Hilhorst, Physica 

A 251, 179 (1998). 
[27] H. Takayasu and A. Yu. Trctyakov, Phys. Rev. Lett. 68, 

3060, (1992). 
[28] I. Jensen, Phys. Rev. E. 50, 3623 (1994). 
[29] P. Grassberger, F. Krause and T. von der Twer, J. Phys. 

ArMath.Gen., 17, L105 (1984). 
[30] N. Menyhard, J. Phys. A 27, 6139 (1994), 
[31] M. H. Kim and H. Park, Phys. Rev. Lett. 73, 2579, 

(1994) . 

[32] K. E. Bassler and D. A. Browne, Phys. Rev. Lett. 77, 
4094 (1996). 

[33] H.Park and H.Park, Physica A 221, 97 (1995). 

[34] N. Menyhard and G. Odor, J. Phys. A 29, 7739 (1996). 

[35] H. Hinrichsen, Phys. Rev. E 55, 219 (1997). 

[36] G. Odor and N. Menyhard, Phys. Rev. E 57,5168 (1998). 

[37] M. Henkel and U. Schollwock, J. Phys. A 34, 3333 (2001). 

[38] K. Park and I. Kim, Phys. Rev E 66, 027106 (2002). 

[39] R. Dickman and M. A. F. de Menenzes, Phys. Rev. E 66, 

045101 (2002). 
[40] G. Odor, cond-mat/0209287. 
[41] P. Grassberger, Z. Phys. B 47, 365 (1982). 
[42] I. Jensen, J. Phys. A 32, 5233 (1999). 



TABLE I. Logarithmic fitting results by the form (20) for 
the one dimensional 3A — » 5A, 2A — > model. 







D 


Pc 


a 


b 


a 


0.2 
0.8 


0.4627 
0.22405 


0.115 
14.12 


3.5 x 10" 5 
-1.001 


0.217 
0.224 
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